#!/bin/sh
# test function definitions for dalton

# check_gen_input:
# generate several input and description files.
# $1 - test name
# $2 - test info.
# $3 - molecule file.
# $4 - dalton file.
# $5 - list of test check routines.
check_gen_input() {
testname=`basename $1`
if [ -x /bin/ksh ]; then
   CHECK_SHELL='#!/bin/ksh'
elif [ -x /bin/bash ]; then
   CHECK_SHELL='#!/bin/bash'
else
   CHECK_SHELL='#!/bin/sh'
fi

# Prefer gawk - we know exactly what it can do.
# awk on Sun does not support functions, need to use nawk for this
if gawk '{print 1}'</dev/null > /dev/null 2>&1; then
   AWK=gawk
elif nawk '{print 1}'</dev/null > /dev/null 2>&1; then
   AWK=nawk
else
   AWK=awk
fi
    if [ "$5" = "" ]; then echo "empty test check file. Test $1 broken."; fi
    echo "$2" > $testname.info
    echo "$3" > $testname.mol
    echo "$4" > $testname.dal
    echo "$CHECK_SHELL" >"$testname.check"
    echo "AWK=$AWK" >> "$testname.check"
    cat >>"$testname.check" <<'%EOF%'
log=$1
ERRLIST=""
. ./functions || { echo "current dir is not the test dir" ; exit 1; }
%EOF%

    echo "$5" >> "$testname.check"
    cat >>"$testname.check" <<'%EOF%'
if [ "$ERRLIST" = "" ]; then
  echo TEST ENDED PROPERLY
  exit 0
else
  echo "Failed tests: $ERRLIST"
  echo THERE IS A PROBLEM
  exit 1
fi
%EOF%
}

# ###################################################################
# miscellaneous tests listed here.

# check_basis:
# check whether number of primitives and contraction parameters
# are as expected. Unlikely to fail but...
# $1 - name of the basis set
# $2 - string like "H 1 1 4s1p 2s1p" containing
#      atom label, number of atoms with this label, charge, 
#      number of primitives, number of contracted functions,
#      types of primitives and contracted functions.
check_basis() {
if [ "$2" = "" ]; then echo "not enough parameters to $0"; return 1; fi
$AWK -v "s=$2" '
/Basis set used is /{
split(s,a)
getline;getline;getline
for(i=1; a[i] != ""; i+=4) {
getline;
if(a[i] != $1 || a[i+1] != $2 || a[i+2] != $3 || a[i+3] != $6) {
    print; print a[i], a[i+1], a[i+2], a[i+3]
    exit(1)
}
}
done=1
}
END {exit(done ? 0 : 1)}' $log || ERRLIST="$ERRLIST; BASIS-SET"
}

# check_scf_energy:
# $1 - calculation type: HF|DFT
# $2 - expected energy
# $3 - max deviation
check_scf_energy() {
[ "$3" = "" ] &&  echo "not enough parameters to check_scf_energy."
$AWK 'function e(v,r){#print "scf_energy: checking ",v, r;
return(v>r||v<-r)}
/Final *'$1' *energy:/{c=1;exit (e($4-('$2'),'$3'))} END{if(!c) exit(1)}
' $log || ERRLIST="$ERRLIST; SCF energy"
}

# check_mp2srdft_energy:
# $1 - calculation type: MP2-srDFT 
# $2 - expected energy
# $3 - max deviation
check_mp2srdft_energy() {
[ "$3" = "" ] &&  echo "not enough parameters to check_mp2_energy."
$AWK 'function e(v,r){#print "mp2_energy: checking ",v, r;
return(v>r||v<-r)}
/= short-range '$1' second order energy         :/{c=1; exit (e($8-('$2'),'$3'))} END{if(!c) exit(1)} 
' $log || ERRLIST="$ERRLIST; MP2srdft energy"
}

# check_dipole:
# $1 - expected dipole moment (au)
# $2 - max deviation
check_dipole() {
[ "$2" = "" ] &&  echo "not enough parameters to check_dipole."
$AWK '/Dipole moment$/{getline;getline;getline;getline;
s=($1<'$1-$2'||$1>'$1+$2'); exit 0;}END{exit (s?1:0)}' $log || ERRLIST="$ERRLIST; Dipole moment"
}

# check_dip_components:
# check components of the first dipole moment encountered in the 
# log file.
# $1 - x component
# $2 - y component
# $3 - z component
# $4 - max deviation
check_dip_components() {
[ "$4" = "" ] &&  echo "not enough parameters to check_dip_components."
$AWK 'function e(a,b,r){return(a-b>r||a-b<-r)}
/Dipole moment components/{
for(i=0;i<5;i++)getline; 
if (s=e($3,'$1,$4')) exit 1;getline;
if (s=e($3,'$2,$4')) exit 1;getline;
if (s=e($3,'$3,$4')) exit 1;exit 0}
END{exit (s?1:0)}' $log || ERRLIST="$ERRLIST; Dipole moment"
}

# check_polarizability:
# check linear response polarizability
# $1 - value of X component
# $2 - value of Y component
# $3 - value of Z component
# $4 - max deviation.
check_polarizabilities() {
[ "$4" = "" ] &&  echo "not enough parameters to check_polarizabilities."
$AWK 'function e(a,b,r){return(a-b>r||a-b<-r)} BEGIN {s=0}
/^@.*XDIPLEN/ {if(!e($8,'$1,$4'))s++}
/^@.*YDIPLEN/ {if(!e($8,'$2,$4'))s++}
/^@.*ZDIPLEN/ {if(!e($8,'$3,$4'))s++}
END{exit (s==3?0:1) }' $log || ERRLIST="$ERRLIST; Polarizabilities"
}

# check_hyperpolarizability:
# check quadratic response/hyperpolarizability
# $1 - value of frequency B
# $2 - value of frequency C
# $3 - combinations of operators, eg. 'X;Y,Y'
# $4 - expected value
# $5 - max deviation.
check_hyperpolarizabilities() {
[ "$5" = "" ] &&  echo "not enough parameters to check_hyperpolarizabilities."
$AWK 'function e(a,b,r){return(a-b>r||a-b<-r)} BEGIN {s=0}
/^@ B-freq/ && $4=='$1' && $7=='$2' && $8=="beta('$3')"{if(!e($10,'$4,$5'))s++}
END{exit (s==1?0:1) }' $log || ERRLIST="$ERRLIST; Hyperpol-ties $4 ($1 $2, $3)"
}

# check_second_hyperpolarizability
# check cubic response/second hyperpolarizability
# $1 - value of frequency B
# $2 - value of frequency C
# $3 - value of frequency D
# $4 - combinations of operators, eg. 'X;Y,Y,X'
# $5 - expected value
# $6 - max deviation
check_second_hyperpolarizabilities() {
[ "$6" = "" ] &&  echo "not enough parameters to check_hyperpolarizabilities."
$AWK 'function e(a,b,r){return(a-b>=r||a-b<=-r)} BEGIN {s=0; isthere=0}
/B-freq:  / { if($3=='$1') {getline;} else {next;}
              if($3=='$2') {getline;} else {next;}
              if($3=='$3') {isthere=1;}
            }
/'$4'/ && isthere {if(!e($3,'$5,$6'))s++}
END{exit (s==1?0:1) }' $log || ERRLIST="$ERRLIST; expected gamma($4)_($1,$2,$3) = $5"
}

# check_general_lr:
# check linear response property P=<<X,Y>>
# $1 - label of X
# $2 - label of Y
# $3 - expected value of P
# $4 - max deviation.
check_general_lr() {
[ "$4" = "" ] &&  echo "not enough parameters to check_general_lr."
$AWK 'function e(a,b,r){return(a-b>r||a-b<-r)} BEGIN {s=0}
/^@ -<< '$1'.*'$2'/ { value=substr($0,31)
                 gsub("D","E",value);
                 if(!e(value,'$3','$4')) s++; 
                 else { print substr($0,31) " " '$3'; }
               }
END{exit (s==1?0:1) }' $log || ERRLIST="$ERRLIST; Linear response"
}

# check_exc_energy:
# check excitation energy
# $1 - symmetry of the excited state
# $2 - number of the excited state
# $3 - excitation energy.
# $4 - max deviation.
check_exc_energy() {
[ "$4" = "" ] &&  echo "not enough parameters to check_exc_energy."
$AWK -v "symm=$1" -v "nr=$2" -v "ene=$3" -v "thr=$4" '
function e(val,r){return(val>r||val<-r)}
/^@ Excited    state    symmetry/{ rsy= ($5==symm)}
/^@ State no:/                   {rst=  ($4==nr)}
rsy && rst && /^@ Excitation energy/{
getline;res=e(ene-$2,thr)
}
END{exit (res) }' $log || ERRLIST="$ERRLIST; Excitation energy"
}

# check_aba_polarizability:
# check the static polarizability as reported by ABACUS property module.
# $1 - the threshold
# $2..$4 -- Ex, Ey, Ez in a.u.
check_aba_polarizability() {
if [ "$4" = "" ]; then
    echo "not enough parameters to $0."; return 1
fi
$AWK -v "thr=$1" -v "ex=$2" -v "ey=$3" -v "ez=$4" '
function e(v,r){return(v>r||v<-r)}
/Static polarizabilities .au/{
for(i=0;i<4;i++) getline;
getline; if($1 != "Ex" || e($2-ex,thr)) {print $2 " " ex;exit(1);}
getline; if($1 != "Ey" || e($2-ey,thr)) {print $2 " " ex;exit(1);}
getline; if($1 != "Ez" || e($2-ez,thr)) {print $2 " " ex;exit(1);}
s=1
}
END{exit(!s)}
' $log || ERRLIST="$ERRLIST; static polarizability"
}

# ###################################################################
# Geometry optimization functions.

# check_mol_gradient:
# check N:th molecular gradient in the file for given atom.
# $1 - which occurence of molecular gradient should be checked.
# $2 - 6-character atom label including symmetry equivalence number (if any)
# $3:$5 - gradient components.
# $6 - tolerance.
check_mol_gradient() {
if [ "$6" = "" ]; then
    echo "not enough parameters to check_mol_gradient."; return 1
fi
$AWK 'function e(v,r){#print "mol_grad: checking " v,",", r;
return(v>r||v<-r)}
 BEGIN {s=0;c=0}
/Molecular gradient/{
if(++c == '$1') { getline; getline;
#skip norms
cemty=0;
while(cempty<9 && substr($0,2,6)!="'"$2"'") {
  if($0=="") cempty++;
  getline;
}
if(cempty==9) { print "mol_gradient: Label '"$2"' not found!"; exit 1;}
split(substr($0,8),a);
exit( e(a[1]-('$3'),'$6') || e(a[2]-('$4'),'$6') || e(a[3]-('$5'),'$6'))
}
}
' $log || ERRLIST="$ERRLIST; molecular gradient"
}

# check_dip_gradient:
# check molecular gradient of dipole moment for given atom
# $1 - number of atoms to check
# $2 - tolerance
# $3 - components following format
# <label1> x <dxx> <dyx> <dzx>
# <label1> y <dxy> <dyy> <dzy>
# <label1> z <dxz> <dyz> <dzz>
# <label2> x <dxx> <dyx> <dzx>
# ...., ie. same as in the output file.
check_dip_gradient() {
if [ "$3" = "" ]; then
    echo "not enough parameters to check_dip_gradient."; return 1
fi
p=`echo $3`
$AWK 'function e(v,r){#print "mol_grad: checking " v,",", r;
return(v>r||v<-r)}
 BEGIN {s=0;c=0}
/Dipole moment gradient \(au\)/ {
getline; # ---------------------------
getline; #
getline; #  Ex             Ey             Ez
getline; #
atoms='$1'; tol='$2';
split("'"$p"'",arr," ");
idx=1;
failed=0;
for(a=0; a<atoms; a++) {
for(cor=0;cor<3; cor++) {
  lbl=arr[idx++]; coor=arr[idx++];
  x=arr[idx++];y=arr[idx++]; z=arr[idx++];
  getline;
  if($1 != lbl || $2 != coor ||
     e(x-$3,tol) || e(y-$4,tol) || e(z-$5,tol)) {
     print $0;
     print "match on ",lbl,":",coor,":", x,":", y,":", z," failed";
     failed++;
  }
}
  getline; #empty
}
}
END {exit(failed)}
' $log || ERRLIST="$ERRLIST; Dipole moment gradient"
}

# check_mol_hessian:
# check the molecular hessian.
# $1 - tolerance
# $2 - the hessian
# checking is tricky because of the large amount of data to be
# verified. We compare all the numbers in the string with given
# tolerance.  When passing string to awk, we compress spaces and
# trailing zeros because some variants of awk (notably on AIX) cannot
# handle strings longer than 399 characters.
check_mol_hessian() {
if [ "$3" = "" ]; then
    echo "not enough parameters to check_mol_hessian."; return 1
fi
p=`echo $3|sed 's,[ 0]* , ,g'`
$AWK 'function e(v,r){#print "mol_hessian: checking " v,",", r;
return(v>r||v<-r)}
 BEGIN {s=0;c=0}
/Molecular Hessian \(au\)/ {
getline; # ---------------------------
getline; #
tol='$1';
split("'"$p"'",arr," ");
idx=1;
failed=0;
while(arr[idx] != "") {
getline; split($0,l," ");
for(i=1; l[i]!= ""; i++) {
 #print "Expected:",arr[idx]," Found:",l[i]
 if(l[i] ~ /-?[.0-9]+/) {
   if(e(arr[idx]-l[i],tol)) failed++;
  } else {
   if(arr[idx] != l[i])     failed++
  } 
  idx++
}
if(failed>5) break
}
}
END {exit(failed)}
' $log || ERRLIST="$ERRLIST; Molecular Hessian"
}

# check_final_geometry:
# $1 - 6-character atom label including symmetry equivalence number
# $2:$4 - coordinates.
# $5    - tolerance
check_final_geometry() {
if [ "$5" = "" ]; then
    echo "not enough parameters to check_final_geometry."; return 1
fi
$AWK 'function e(v,r){#print "final_geometry;checking ",v, r;
return(v>r||v<-r)}
BEGIN {s=0;c=0}
/Final geometry/{ getline; getline;getline;
while($0 != "" && substr($0,2,6)!="'"$1"'") getline;
if($0 == "") { print "final_geometry: Label '"$1"' not found!"; exit 1;}
split(substr($0,8),a);
exit( e(a[1]-('$2'),'$5') || e(a[2]-('$3'),'$5') || e(a[3]-('$4'),'$5'))
}
' $log || ERRLIST="$ERRLIST; final geometry"
}

# check_final_energy:
# $1 - final geometry optimization energy
# $2 - tolerance.
check_final_energy() {
if [ "$2" = "" ]; then
    echo "not enough parameters to check_final_energy."; return 1
fi
$AWK 'function e(v,r){#print "checking ",v, r;
return(v>r||v<-r)}
/Energy at final geometry is/ { exit(e($7-('$1'),'$2'))  }
' $log || ERRLIST="$ERRLIST; final energy"
}

# check_rot_constants:
# check rotational constant
# $1 - expected value in MHz
# $2 - expected value in cm-1 
# Tested for a linear molecule.
check_rotational_const() {
if [ "$2" = "" ]; then
    echo "not enough parameters to check_rotational_const."; return 1
fi
$AWK '
/ B = /{if($3!='$1'||$6 != '$2')failed++}
END    {exit(failed)}
'  $log || ERRLIST="$ERRLIST; ROTATIONAL CONSTANTS $1"
}

# check_vibr_freqs:
# check the vibrational frequencies and intensities.
# $1 - tolerance
# $2 - the IR freq and intensity data
# use the same method as in molecular hessian. 
check_vibr_freqs() {
if [ "$2" = "" ]; then
    echo "not enough parameters to check_vibr_freq."; return 1
fi
p=`echo $2`
$AWK -v "tol=$1" 'function e(v,r){#print "vibr_freq: checking " v,",", r;
return(v>r||v<-r)}
 BEGIN {s=0;c=0}
/ Vibrational Frequencies and IR Intensities/ {
for(i=0; i<9; i++) getline;
split("'"$p"'",arr," ");
idx=1;
failed=0;
for(idx=1;arr[idx] != ""; idx+=6) {
getline; split($0,l," ");
for(i=1; i<=6; i++) {
 #print "Expected:",arr[idx+i-1]," Found:",l[i]
 if(l[i] ~ /-?[.0-9]+/) {
   if(e(arr[idx+i-1]-l[i],tol)) {failed++;print "expected",arr[idx+i-1]," Found ", l[i],tol}
  } else {
   if(arr[idx+i-1] != l[i])     {failed++;print "b"}
  } 
}
if(failed>5) break
}
}
END {exit(failed)}
' $log || ERRLIST="$ERRLIST; Vibrational Frequencies"
}

# check_exp_grad:
# Check the orbital gradient.
# $1 - number of the basis function
# $2 - exponent of the basis function (for input cross-checking)
# $3 - expected value of the gradient
# $4 - tolerance
check_exp_grad() {
if [ "$4" = "" ]; then
    echo "not enough parameters to check_exp_grad."; return 1
fi
$AWK 'function e(v,r){#print "checking ",v, r;
return(v>r||v<-r)}
$1== "##" && $2 == '$1' && $7=='$3' { exit(e($8-('$3'),'$4'))  }
' $log || ERRLIST="$ERRLIST; ORB. GRADIENT $1"
}

check_magnetizabilities() {
return 0
}

# check_mag_tensors:
# check magnetizability tensors, total, diamagnetic and paramagnetic
# contributions.
# $1 - threshold
# $2 - total tensor
# $3 - nuclear contribution to the tensor
# $4 - electronic contribution to the tensor
check_mag_tensors() {
if [ "$4" = "" ]; then
    echo "not enough parameters to check_mag_tensors."; return 1
fi

$AWK -v "tol=$1" -v "t=`echo $2`" -v "d=`echo $3`" -v "p=`echo $4`" '
function e(v,r){return(v>r||v<-r)}
/Total magnetizability tensor/ {
split(t,a); for(i=0;i<5;i++)getline;
if(e($2-a[1],tol) || e($3-a[2],tol) ||e($4-a[3],tol)) exit 1
getline;
if(e($2-a[4],tol) || e($3-a[5],tol) ||e($4-a[6],tol)) exit 1
getline;
if(e($2-a[7],tol) || e($3-a[8],tol) ||e($4-a[9],tol)) exit 1
foundt=1
}
/Diamagnetic magnetizability tensor/ {
split(d,a); for(i=0;i<5;i++)getline;
if(e($2-a[1],tol) || e($3-a[2],tol) ||e($4-a[3],tol)) exit 1
getline;
if(e($2-a[4],tol) || e($3-a[5],tol) ||e($4-a[6],tol)) exit 1
getline;
if(e($2-a[7],tol) || e($3-a[8],tol) ||e($4-a[9],tol)) exit 1
foundd=1
}
/Paramagnetic magnetizability tensor/ {
split(p,a); for(i=0;i<5;i++)getline;
if(e($2-a[1],tol) || e($3-a[2],tol) ||e($4-a[3],tol)) exit 1
getline;
if(e($2-a[4],tol) || e($3-a[5],tol) ||e($4-a[6],tol)) exit 1
getline;
if(e($2-a[7],tol) || e($3-a[8],tol) ||e($4-a[9],tol)) exit 1
foundp=1
}
END{exit( (foundt && foundd && foundp) ? 0 : 1) }
' $log  || ERRLIST="$ERRLIST; MAGNETIZABILITY TENSORS"
}

# check_g_tensors:
# check rotational g-tensors, total, nuclear and electronic
# contributions.
# $1 - threshold
# $2 - total tensor
# $3 - nuclear contribution to the tensor
# $4 - electronic contribution to the tensor
check_g_tensors() {
if [ "$4" = "" ]; then
    echo "not enough parameters to check_g_tensors."; return 1
fi

$AWK -v "tol=$1" -v "t=`echo $2`" -v "nu=`echo $3`" -v "el=`echo $4`" '
function e(v,r){return(v>r||v<-r)}
/Molecular rotational g-tensor in principal axis system/ {
split(t,a); for(i=0;i<3;i++)getline;
if(e($2-a[1],tol) || e($3-a[2],tol) ||e($4-a[3],tol)) exit 1
getline;
if(e($2-a[4],tol) || e($3-a[5],tol) ||e($4-a[6],tol)) exit 1
getline;
if(e($2-a[7],tol) || e($3-a[8],tol) ||e($4-a[9],tol)) exit 1
foundt=1
}
/Nuclear contribution in principal axis system/ {
split(nu,a); for(i=0;i<3;i++)getline;
if(e($2-a[1],tol) || e($3-a[2],tol) ||e($4-a[3],tol)) exit 1
getline;
if(e($2-a[4],tol) || e($3-a[5],tol) ||e($4-a[6],tol)) exit 1
getline;
if(e($2-a[7],tol) || e($3-a[8],tol) ||e($4-a[9],tol)) exit 1
foundn=1
}
/Electronic contribution in principal axis system/ {
split(el,a); for(i=0;i<3;i++)getline;
if(e($2-a[1],tol) || e($3-a[2],tol) ||e($4-a[3],tol)) exit 1
getline;
if(e($2-a[4],tol) || e($3-a[5],tol) ||e($4-a[6],tol)) exit 1
getline;
if(e($2-a[7],tol) || e($3-a[8],tol) ||e($4-a[9],tol)) exit 1
founde=1
}
END{exit( (foundt && foundn && founde) ? 0 : 1) }
' $log  || ERRLIST="$ERRLIST; G-TENSORS"
}
